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Abstract 

The Born effective charges of component atoms and phonon spectra of 
a tetrahedrally coordinated crystalline ice are calculated from the first 
principles method based on density functional theory within the gen- 
eralized gradient approximation with the projected augmented wave 
method. Phonon dispersion relations in a 3 x 1 x 1 supercell were eval- 
uated from Hellmann-Feynman forces with the direct method. This 
calculation is an additional work to the direct method in calculating 
the phonon spectra which does not take into account the polariza- 
tion charges arising from dipole interaction of molecules of water in 
ice. The calculated Born effective polarization charges from linear re- 
sponse theory are supplied as the correction terms to the dynamical 
matrix in order to further investigate the LO-TO splitting of the po- 
lar modes of ice crystal at k = which has long been speculated for 
this system especially in the region between 28 and 37 meV both in 
the theoretical and experimental studies. Our results clearly show the 
evidence of splitting of longitudinal and transverse optic modes at the 
k = 0-point in agreement with some experimental findings. 

Keywords: Ice, Density functional calculations, Phonon spectra, Born 
effective charges 

1. INTRODUCTION 

The vibrational studies of solid ice have gained attraction for scientific in- 
vestigation over the decades. Experimentally, the vibrational spectra of the 
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different phases of ice have been investigated using infra-red (IR) and Ra- 
man techniques, but due to the proton disordering in most such structures, 
the normal selection rules governing the interaction of radiation with these 
lattices are broken and hence any analysis of the spectra is difficult. On 
the other hand, the IR and Raman spectra are very sensitive to the in- 
tramolecular modes involving O-H stretching and bending, and less sensitive 
to the intermolecular modes involving the vibrations where the whole water 
molecules are moving against each other. Therefore, in normal circumstances, 
only limited information (the acoustic frequencies in particular) can be ob- 
tained in the translational region. Inelastic incoherent neutron scattering 
(HNS) is only more reliable because its spectrum is directly proportional to 
the phonon density of states weighed by mean square amplitude associated 
with each mode and also the selection rule is not involved as all modes are 
measured simultaneously Pfl |2] . Despite the detail information that can be 
obtained from the present days methods of coherent and incoherent scatter- 
ing on the vibrational motions of the atoms or molecules of most systems 
based on the peak positions, intensities and their width, there are still many 
problems about the ice systems. These problems as mentioned, are usually 
associated with the the proton disordering. The phonon spectra measure- 
ments of Renker j3j on Ice-Ih (D2O) made by time of flight using a chopper 
spectrometer were extensive but not complete in that information above 20 
meV in the (001) direction and above 30 meV in other directions is missing, 
presumably because of a lack of scattering intensity and of the nature of the 
disordering of the protons in the ice structure. 

First-principles calculation of crystal structure and crystal properties is 
becoming standard technique, and the progress in the methods, algorithms, 
and computer capabilities allows to study larger systems of solids crystals 
in which crystalline ice cannot be left out. The method has recently gained 
ground not only because of its reliability in the study of static and dynamics 
properties of ice |U |3] but also some important features such as modeling 
ordered periodic ice structure |E] , and also to probe the nature of hydrogen 
bond in different geometries £[] . Its theoretical counterpart such as the clas- 
sical modeled potential through empirical method [HI E] had some success 
in describing some important dynamical features, but to date, there is none 
capable of describing the ice dynamics and related properties across its whole 
spectral range and describing certain key spectral features. 

There are two techniques currently in use in the first-principles method 
in the study of lattice dynamics of crystals: The linear response method and 
the direct method. In the linear response method the dynamical matrix is 
obtained from the modification of the electron density, via the inverse di- 
electric matrix. The dielectric matrix is calculated from the eigenfunctions 
and energy levels of the unperturbed system jTU|- R can be determined at 
any wave vector in the Brillouin zone with the computational effort required 
comparable to that of a ground state optimization. Only linear effects, such 
as harmonic phonons, are accessible to this technique. On the other hand, 
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the direct-method is based on the solution of the Kohn-Sham equation and 
it allows one to study both linear and non-linear effects. The calculations 
deal with a supercell, which allows explicit account to be taken of any per- 
turbation. This method is rather straightforward computationally and there 
are a few standard software packages. Within the direct method the phonon 
frequencies are calculated from Hellmann-Feyman forces generated by the 
small atomic displacements, one at a time. Hence using the information 
of the crystal symmetry space group the force constants are derived, and 
the dynamical matrix is built and diagonalized, and its eigenvalues arranged 
into phonon dispersion relations. In this way, phonon frequencies at selected 
high frequencies at high-symmetry points of the Brillouin zone can be cal- 
culated [TT]. However, when the interaction range ceases to be within the 
supercell, phonons at all wave vectors are determined exactly. The above 
statement has to be modified for polar crystals for which the macroscopic 
electric field splits off the infrared-active optic modes. The long-range part 
of the Coulomb interaction corresponds to the macroscopic electric field aris- 
ing from ionic displacements. 

Ice is a tetrahedrally covalently bonded polar system whose dipole-dipole 
interactions give rise to the electric field when they are disturbed. The origin 
of the splitting is therefore the electrostatic field created by long wavelength 
modes of vibrations in such crystals. Usually a microscopic electric field in- 
fluences only the LO modes while TO modes remain unaltered. The field 
therefore breaks the Born-von Karman conditions, as a consequence with a 
direct method only finite wave vector k ^ calculations are possible. The 
LO/TO splitting has therefore been found by calculating the effective Born 
charge tensor and electronic dielectric constant introduced into the dynam- 
ical matrix in the form of a non-analytical term ^2] or by calculating LO 
modes from elongated supercells jTSI as will be discussed below. 

In our previous work, we have applied with success the direct method to 
the calculation of phonon dispersion of ice and its corresponding vibrational 
density of states [Hj. The method reproduces the important features in the 
translational mode, librational mode, bending as well as the stretching re- 
gion in comparison to the experimental results. The only ingredient that 
was missing in the previous work has been explained above, i.e., a macro- 
scopic electric field arising from dipolar interaction which is not taken into 
account. Despite the huge computational demand of this problem, we did 
the additional calculation of Born-effective charges which is supplied as in- 
gredients to the direct method to identify the long-range part of interatomic 
force constants and makes the interpolation of phonon frequencies tractable. 
Our overall aim is to help resolve the discrepancies in the reported phonon 
frequencies especially the puzzle behind the LO/TO splitting of some opti- 
cal modes at k = and provide first principles Born-effective charges and 
dielectric tensors for direct method phonon calculations. 
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2. METHOD OF CALCULATION 



The calculation of phonon dispersion relations were performed with the di- 
rect method. The direct method uses the Hellmann-Feymann (HF) forces 
calculated for the optimized supercell with one atom displaced from equilib- 
rium position, derived from the force constants using the symmetry elements 
of the space group of the crystal, and calculates phonon frequencies by diag- 
onalizing the dynamical matrix. 

In this work, the ice crystal structure optimization and calculation of HF 
forces have been performed with the density functional theory using the 
PAW method within the generalized gradient approximation (GGA), as im- 
plemented in the Vienna Ab Initio Simulation Package (VASP) jT^j software. 
A unit cell of ice crystal was prepared in a cubic box according to Fig.^with 
8 molecules of water. All the atomic degrees of freedom were relaxed with 
high precision. The optimum Monkhorst Pack of 4 x 4 x 4 fc-point was used in 
addition to the GGA of Perdew-Wang to describe the exchange-correlation 
and the hydrogen bonding of water. We used a energy cut-off of 500 eV 
because the 2p valence electrons in oxygen require a large plane wave basis 
set to span the high energy states described by the wavefunction close to 
the oxygen nucleus, and also the hydrogen atoms require a larger number of 
planes waves in order to describe localization of their charges in real space. 
The relaxed geometry for the unit cell from the initial configurations contain- 
ing 8 molecules is shown in Fig. ^ The starting geometry of the molecules in 
the cubic simulation box shown is such that no hydrogen bonds were present 
but the positions of oxygen atoms follow the tetrahedral geometry. After the 
relaxation, all the protons perfectly point to the right direction of oxygen 
atoms and make the required hydrogen bonds necessary as indicated by the 
dotted lines in Fig. ^ to preserve the tetrahedral orientation of the ice struc- 
ture. This final structure is in accordance to Bernal Fowler's rules [HI] which 
are based on the statistical model of ice. 

The relaxed geometry is tetragonal with calculated lattice parameters a = 
6.1568 A, b = 6.1565 A, c = 6.0816A i.e. with c/a ratio « 0.988. The ex- 
perimental lattice constant reported by Blackman et. al. ^7] is 6.35012652 
A for the cubic geometry. The calculations of force constants was carried 
out by considering a 3 x 1 x 1 supercell containing 24 molecules of water 
which is obtained by matching 3 tetragonal unit cells. At the first step of the 
calculation, the PHONON software is used to define the appropriate crystal 
supercell for use of the direct method. The phonon frequencies cj(k,j) are 
calculated as square roots of eigenvalues of the supercell dynamical matrix: 



where the e(k, j) are the polarization vectors. The supercell dynamical ma- 
trix is defined as 
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x exp(— 27rik.[R(0, fi) — R(m, u)]), 

where the summation over m runs over all atoms of the supercell; M M , M v 
and R(0, /i), R(m, v) are atomic masses and equilibrium vectors, respectively; 
and k is the wavevector. The cummulant force constants <f>fp are the sums 
of terms containing the second derivatives of the ground-state energy with 
respect to the position vectors of interacting atoms i and j. The HF forces 
in the direct method are derived using 

Fi(n,u) = - ^ C (^ z/ ; m ) / i ) M i( m ) / i )> ( 3 ) 

m,i/,j 

where Uj(m, fi) is an amplitude of displacement of an atom in the supercell 
specially shifted from the equilibrium position. 

The symmetry of the supercell and the site symmetry of the non-equivalent 
atoms usually considerably reduce the number of displacements needed for 
reconstruction of 3>^ c - We are unfortunate in our case because of the hy- 
drogen bonding fluctuations which makes the crystal structure complicated. 
Therefore, we have to consider the whole 24 atoms in the supercell as inde- 
pendent displacements: In the positive and negative non-coplanar, x, y and 
z directions. As done for the primitive unit cell, all the internal coordinates 
were relaxed until the atomic forces were less than 10 -4 eV/A. Complete in- 
formation of the values of force constants were obtained by displacing every 
atom of the primitive unit cells by 0.02 A in both positive and negative x, y 
and z directions. Therefore, minimization of the anharmonic effects and sys- 
tematic errors are achieved by calculating with Eq. fusing forces arising 
from both positive and negative displacements Uj. As mentioned above, we 
use a 3 x 1 x 1 supercell, which implies that 3 points in the direction [100] are 
treated exactly according to the direct method. The points are [£00] , with ( 
= 1, 1/3, 2/3. We calculate forces induced on all atoms of the supercell when 
a single atom is displaced from its equilibrium position, to obtain the force 
constant matrix, and hence the dynamical matrix. This is then followed by 
diagonalization of the dynamical matrix which leads to a set of eigenvalues 
for the phonon frequencies and the corresponding normal-mode eigenvectors. 
The vibrational density of states (VDOS) is obtained by integrating over k- 
dependent phonon frequencies from the force- const ant matrix in supercells 
derived from the primitive molecule unit cells. 

For the ionic crystals the macroscopic electric field is taken into account by 
adding to Eq. (Q) the non-analytic term of the dynamical matrix at the wave 
vector k = |18j . However, since one knows the phonon frequencies only 
at discrete wave vectors, it is justified to extend the non-analytical term to 
the k ^ region, through multiplying it by the Gaussian damping factor. 
Therefore we replace Eq. (Q) by the following expression: 

D^(k;H = D sc ai/3 (k;H 

47re 2 [k-Z*(//)] a [k-Z*(^ 
Ve^M^Mv |k| 2 
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x exp[— 27rig • (r(/x) — r(u))} 
[/u \ 2 

x d(q)exp I —n 2 




(4) 



where k is the wave vector within the Brillouin zone with its centre at the 
reciprocal-lattice vector g, V stands for the volume of the primitive unit 
cell, and M M , r M are atomic masses and internal positions, respectively. The 
Z* (/i) are the tensors of the Born-effective charges, eoo is the electronic part 
of the dielectric constant and p (x, y and z) are damping factors; then the 
non-analytical term vanishes close to the zone boundary. Consideration of 
the effective charges leads to the LO/TO splitting of the optical parts of the 
phonon modes of of ice at the T-point as discussed in the next Section. This 
observation has long been been speculated both from theory and experiment. 
By definition ^3], the Born effective charge tensor Z*j iQ:/ g quantifies to linear 
order the polarization per unit cell (P) generated by zone-center k=0, created 
along the direction f3 when the atoms of sublattice i are displaced in the 
direction a under the condition of zero electric field. It is calculated according 
to the equation: 



i,af3 



dPg 



The macroscopic dielectric constant is found via the relation 



1 + 



4ttP 



E 



(5) 



(6) 



where E = E ext — 4nP is the total macroscopic electric field. 



3. RESULTS AND DISCUSSION 

Table I shows the dielectric constants of the tensor and the Born-effective 
charges calculated according to Eq. (JHJ) and (JSJ) implemented in PWSCF 2 .i 
code j20j. The dielectric constant tensor is symmetric with non-zero off- 
diagonal terms but with negligibly contribution in comparison to the diago- 
nal element. The diagonal term ~ 1.88 is in a very good range for the high 
frequency limit (Thz) of the dielectric constant of ice [21] . Under an applied 
field the individual molecules are polarized by the field. This involves the 
displacements of the electrons relative to the nuclei and small distortions of 
the molecules under the restoring forces. The response to a change in field is 
very rapid, so that the effects are independent of frequency up to microwave 
frequencies. The polarization in ice in general is due to the reorientation of 
molecules or bonds, that is, the energies of some of the proton configurations, 
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that are compatible with the ice rules [TT)|l22j. are lowered relative to others, 
so that in thermal equilibrium there is net polarization of ice. The achieve- 
ment of this equilibrium state is a comparatively slow process that requires 
thermal activation and local violations of the ice rules. 
Our results for the dynamical effective charges are shown in Table I for 24 
non-equivalent atoms. The complication of the system due to the protons 
re-orientation and the hydrogen bonding make the symmetry consideration 
difficult. The charge neutrality condition [TH] requires that the acoustic mode 
frequencies vanish for k =0 such that 

K# = o- (7) 

Therefore, this condition is satisfied to at least order of 10 -4 electron which 
is accurate enough for any reliable calculation. The unequal effective charge 
tensor components and Zq of each hydrogen and oxygen atoms Z xx ^ 
Z* y 7^ Z* z are due to the broken symmetry arising from the lattice distor- 
tion. The values alternate among the component atoms of H and also for 
O in order to preserve the overall neutrality. For instance, the Z* zz of the 
hydrogen atoms is 0.624(±0.001) electron, while the other two elements Z* y 
and Z* zz alternate within 0.666(±10 -4 ) electron. The off-diagonal elements 
have the values which ranges within ±0.408, ±0.407, ±0.402, ±0.377 and 
±0.382 with deviation ±10 -4 electron. Similar features are observed for the 
oxygen atoms with Z* zz -1.070 ±10~ 4 electrons while Z xx and Z* y alternate 
between -1.084 and -0.961 with deviation (±10~ 4 ) electrons. Some of the 
off-diagonal contributions are too small. The observed anisotropic features 
can be attributed to the complexity of the hydrogen bond during the electron 
transfer process and also due to the dipole interaction of the water molecules 
in Eq. |U 

Figure El shows the dispersion relation obtained by supplying the calculated 
effective charges and the corresponding highest frequency dielectric con- 
stants, shown in Table I, as the correction from the analytical term which was 
added to the dynamical matrix as explained in Section 2. We also compare 
the dispersion obtained in the absence of these charges to see the magnitude 
of splitting in the optical mode. The calculated VDOS for both disper- 
sions do not appreciably change much because the states are not complete 
as it requires summation over all points in the first Brillouin zone. The LO 
modes that was formally degenerate in the absence of Z* at 27.1 meV in the 
translational region is now shifted to a higher value 30.2 meV as shown in 
Fig. El Except for the isotopic effect due to the different masses of hydrogen 
and deuterium atoms in ice (see the experimental dispersion on the right of 
Fig. El). This observation correctly shows the splitting in the optical modes 
due to the dipole interaction of the water molecules which correspondingly 
induce a dipole moment in the optical mode. Apart from the translational 
region, as shown in Fig. 01 other splitting of LO modes occur at 111.0 meV in 
the librational region with a small shift to 112.5 meV. Also, there is a small 
shift of 0.1 meV from 204.6 meV in the bending region. The large shift of 
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about 12.0 meV is observed from 362 meV in the stretching region because 
the strength of dipole interaction is large when there is symmetric stretching 
iyx) of O-H bond, while a tiny effect observed in the antisymmetric stretch- 
ing is due to the compensation arising from simultaneous bond lengthening 
and shortening of the O-H covalent bond. 

4. SUMMARY 

In summary, we have performed the analysis of the lattice dynamics in crys- 
talline tetrahedrally coordinated ice and found a very strong influence of the 
dipole interaction on the phonon spectra in the optical regions of ice. Ab 
initio calculations clearly show the splitting in the region between 27.1 and 
30.2 meV of the translational region . This observation can be correlated 
with the experimental observation using high resolution inelastic neutron 
measurements of the phonon density of states, in which two separate molec- 
ular optical bands at about 28 and 37 meV for ice Ih and ice Ic have been 
observed [2]. Other splitting of LO modes in the librational, bending and 
stretching region were also predicted due to the dipole interactions. The 
calculation also shows the extent to which the direct method can be used 
to calculate the phonon spectra of the dipole system. To our knowledge, a 
strong influence of the dipole interaction on the lattice dynamics of ice from 
effective charges calculation was not yet reported. 
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TABLES 



Table 1 Calculated dielectric constants and Born effective charges of cubic 
ice using linear response in PWSCF 2 .i code |2Uj . 
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FIGURE CAPTIONS 



Fig. ^ Initial and the relaxed geometry of the unit cell of ice. The ice 
structure was initially packed in a cubic unit cell with initial lat- 
tice constant taken from the literature [H] to be 6.35 A. There 
are no hydrogen bonds in the initial prepared structure shown on 
the left but were perfectly formed after the relaxation according 
to the GGA calculation. The relaxed geometry has the values of 
a ~ b 7^ c which implies that the relaxed structure is tetragonal 
with c/a ratio « 0.988 

Fig. [2] Calculated dispersion curves for ice (a) with no effective charges 
Z* taken into account and (b) with Z* taken into account. The 
curve on the left is the integrated phonon density of states. The 
frequencies v\, v^-, are v% are respectively bending, symmetric and 
anti-symmetric stretching analogous to the vibrational mode of 
an isolated water molecule [II) . Note that the vibrational phonon 
density of states is not complete as it requires summation over all 
points in the first Brillouin zone, nevertheless the calculated G(o;) 
for both with and without Z* do not differ. 

Fig. 03 Phonon dispersion in the transitional region showing the split- 
ting of LO mode in the translational region. The dispersions are 
compared for the cases of both with and without the dipole inter- 
action through the calculation of effective charges Z*. The case 
with Z* is marked by a*. The experimental dispersion taken 
from [2]) is shown on the right 

Fig. 21 Phonon dispersion in the librational, bending and stretching re- 
gion showing the splitting of LO modes. The case with Z* is 
marked by a* as in Fig. El 
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Table I: Adeagbo et. al 

Dielectric constant in cartesian axis 

( 1.881163770 -0.000048586 -0.000019997 ) 
( -0.000048586 1.881154759 0.000112163 ) 

( -0.000019996 0.000112161 1.883382267 ) 



Effective charges E-U in cartesian axis 



Water 


molecule (1) 


(hydrogen 


atom 


1) 






0.66721 


0.37702 







40166 




0.40921 


. 53464 







38200 




0.40785 


0.36304 







62434 


Water 


molecule (1) 


(hydrogen 


atom 2) 






0.66611 


-0.37720 







40208 




-0.40909 


. 53493 




-o 


38272 




0.40784 


-0.36353 







62525 


Water 


molecule (2) 


(hydrogen 


atom 


1) 






0.66619 


0.37729 




-o 


40199 




0.40909 


0.53496 




-o 


38262 




-0.40785 


-0.36360 







62507 


Water 


molecule (2) 


(hydrogen 


atom 


2) 






0.66694 


-0.37691 




-0 


40154 




-0.40907 


0.53471 







38197 




-0.40768 


0.36296 







62440 


Water 


molecule (3) 


(hydrogen 


atom 


1) 






0.53461 


0.40937 




-o 


38251 




0.37708 


0.66692 




-o 


40207 




-0.36332 


-0.40814 







62483 


Water 


molecule (3) 


(hydrogen 


atom 2) 






0.53459 


-0.40938 







38266 




-0.37714 


0.66694 




-o 


40229 




0.36338 


-0.40822 







62507 


Water 


molecule (4) 


(hydrogen 


atom 


1) 






0.53516 


-0.40898 




-o 


38225 




-0.37710 


0.66634 







40159 




-0.36333 


0.40759 







62478 


Water 


molecule (4) 


(hydrogen 


atom 


2) 






0.53515 


0.40892 







38210 




0.37715 


0.66626 







40141 




0.36335 


0.40747 







62434 


Water 


molecule (5) 


(hydrogen 


atom 


1) 






0.66666 


-0.37682 




-o 


40144 




-0.40902 


0.53480 







38199 




-0.40761 


0.36297 







62446 
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Water 


molecule (5) 


(hydrogen 


atom 2) 






0.66648 


. 37748 




-0 


40192 




0.40929 


0.53519 




-0 


38256 




-0.40784 


-0.36361 







62469 


Water 


molecule (6) 


(hydrogen 


atom 1) 






0.66694 


0.37695 







40158 




0.40919 


0.53476 







38206 




0.40781 


0.36307 







62443 


Water 


molecule (6) 


(hydrogen 


atom 2) 






0.66639 


-0.37738 







40199 




-0.40926 


0.53512 




-0 


38263 




0.40780 


-0.36352 







62484 


Water 


molecule (7) 


(hydrogen 


atom 1) 






0.53456 


0.40948 




-0 


38257 




0.37720 


. 66723 




-0 


40225 




-0.36340 


-0.40832 







62487 


Water 


molecule (7) 


(hydrogen 


atom 2) 






. 53464 


-0.40928 







38261 




-0.37704 


0.66666 




-0 


40213 




0.36332 


-0.40807 







62503 


Water 


molecule (8) 


(hydrogen 


atom 1) 






0.53514 


0.40901 







38225 




0.37720 


0.66637 







40161 




0.36345 


0.40765 







62469 


Water 


molecule (8) 


(hydrogen 


atom 2) 






0.53515 


-0.40887 




-o 


38209 




-0.37703 


0.66620 
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Figure 1: Adeagbo et al. 
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Figure 2: Adeagbo et al. 
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Figure 3: Adeagbo et al. 
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Figure 4: Adeagbo et al. 
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